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Summary 

The numerical analysis of the incompressible Navier-Stokes 
equations are becoming important tools in the understanding of 
some fluid flow problems which are encountered in research as well 
as in industry. With the advent of the supercomputers, more realistic 
problems can be studied with a wider choice of numerical algorithms. 

This report presents an alternative formulation for viscous 
incompressible flows . The incompressible Navier-Stokes equations 
are cast in a velocity /vorticity formulation (1). This formulation 
consists of solving the Poisson equations for the velocity components 
and the vorticity transport equation. 

Two numerical algorithms for the steady two-dimensional 
laminar flows are presented. The first method is based on the 
actual partial differential equations. This uses a finite-difference 
approximation of the governing equations on a staggered grid. 

The second method uses a finite element discretization with the 
vorticity transport equation approximated using a Galerkin 
approximation and the Poisson equations are obtained using a least 
squares method. The equations are solved efficiently using Newton’s 
method and a banded direct matrix solver (LINPACK). 

The method is extended to steady three-dimensional laminar 
flows and applied to a cubic driven cavity using finite difference 
schemes and a staggered grid arrangement on a Cartesian mesh. The 
equations are solved iteratively using a plane zebra relaxation 
scheme. 

Currently, a two-dimensional, unsteady algorithm is being 
developed using a generalized coordinate system. The equations are 
discretized using a finite-volume approach. This work will then be 
extended to three-dimensional flows. 



Algorithm Development 
Background 


Some of the widely used solution methods in the study of 
viscous incompressible flows are the primitive variable formulation 
for both 2-D and 3-D and the vorticity /stream function for 2-D and 
plane flows. Under the class of the primitive variable formulations 
are the artificial compressiblity method as developed by Kwak , et. al. 
(2) and the fractional step method by Kim, and Moin (3). More 
recently, Rosenfeld, et. al extended the fractional step method for 
unsteady 3-D flow problems (4). One of the major concerns about 
this formulation is the prescription of pressure on the boundary and 
the method of obtaining the pressure. Mass conservation is an 
important criterion that must also be satisfied. While the pressure 
does not explicitly appear in the continuity equation it acts as an 
important parameter in ensuring that continuity is satisfied. 

Pressure can be eliminated from the equation by introducing the 
vorticity/stream function formulation for 2-D . This formulation 
conserves mass automatically. This method can be extended to 3-D 
using the vorticity/vector potential method. However, this extension 
is not straightforward. Vorticity can be eliminated from this 
equation to obtain a single biharmonic stream function formulation. 
One of the issues of this formulation is that one may lose accuracy if 
the formulation is not done correctly. Convergence with this method 
is also slow. 


The velocity/vorticity method is a combination of the primitive 
variable formulation and the stream function formulation. This 
method eliminates the problem with the pressure but introduces the 
problem of the vorticity boundary conditions. Many different 
techniques to overcome the problem of the vorticity boundary 
conditions have been presented in the literature. They are grouped 
into 3 different types of schemes. The first scheme uses the 
relationship between velocity and stream function. (12, 15). 

The second class are those proposed by Quartapelle and Quartapelle 
and Valz-Gris( 16,17). These papers showed that in order for the 
boundary conditions on the velocity be satisfied the vorticity should 
evolve subject to an integral constraint. The third class of method 
uses the vortex blob methods as introduced by Chorin. The present 
treatment of vorticity boundary condition fall under the second class. 


The boundary conditions on the velocity induce a constraint on the 
vorticity. The equation for the definition of vorticity is then solved 
as part of the flow. 


The equations that are being used in this type of approach is of a 
higher order than the original partial differential equations. Because 
of this, less restrictive boundary conditions can be used. This may be 
advantageous for those applications where restrictive boundary 
conditions are problematic. 

Another advantage of this method when used with a staggered mes 
arrangement is that it gives a compact and accurate representation of 
the conservation of mass and vorticity. The natural decoupling of 
the governing equations is also beneficial since this will lead to 
simple treatment of the boundary conditions. This method seems 
favorable for most 2-D flow applications. However, the extension to 
3-D may present some drawbacks. 

One of the major concern is the increased number of equations and 

unknowns. You now have three equations for the velocity and three 

vorticity components. This problem may be alleviated by not 

solving the equations directly as was done in the 2-D case but rather 

by using an efficient relaxation which will be competitive with the 
other formulations. This problem will be addressed later on when 
the 3-D unsteady formulation using a generalized coordinate system 
is finished. 



Governing Equations 

The continuity equation, the definition of vorticity and the 
vorticity transport equation are given by equations (l)-(3). 

v*q = s , v x q = w 

— t o 

= Z • vq + vvZ 

In equation (1), s represents a source distribution in the field. If s 
# 0, then continuity equation cannot be automatically satisfied by 
introducing a stream function for two dimensional and axisymmetric 
flows. Moreover, the use of stream functions in multiply connected 
domains requires special treatment of boundary conditions. 

To solve equations (l)-(2), with prescribed normal velocity 
components on the boundary, Rose (4,5) introduced a compact finite 
difference scheme which leads to an ill-conditioned system of 
equations. Osswald, Ghia and Ghia (6) used a direct solver for their 
discrete velocity equations decoupled from the vorticity equation. 
Recently Chang and Gunzburger(7) proposed a finite element 
discretization of equations (l)-(2), which leads to a rectangular 
algebraic system of equations. 

The solution of equations (1),(2), with prescribed normal velocity 
component on the boundary,exists only if some compatibility 
relations hold. For example, using Gauss theorem one obtains 

JjJ v-q dV = JJJ s dV = IJ q*n dA 

V V A 

where n is the unit vector perpendicular to A. Similarly, it can be 
shown that 

HI vxq dV =* HI w dV * II nxq dA 

V V A 

Also, the following vector identity must hold 

y.vxq = a 0 

Fasel proposed to take the curl of equation (2) 



Substituting equation (1) in equation (7) gives 


-v^q = vxw - vs 


Assuming w is known, the Poisson equations with Dirichlet 
boundary conditions can be solved for the velocity components. In 
general, equation (8) admits more solutions than those of equations 
(l)-(2). For the case of "inviscid" flow, in a simply connected domain, 
Fix and Rose (8) and Phillips (9) showed that the solution of 
equations (1) and (2) with prescribed normal velocity components 
can be obtained as a solution of equation (8). Besides the normal 
velocity components, equation (2) is used as boundary conditions for 
equation (8). With this choice, equation (8) has a unique solution, 
and spurious solutions are excluded. 

To extend the above analysis to the viscous flow case, the 
prescribed velocity at the boundary is decomposed into two sets: the 
normal and the tangential components. The normal component 
together with the definition of the vorticity are imposed as boundary 
conditions for equation (8). The tangential velocity components 
provide boundary conditions for the vorticity transport equation. If 
equation (8) and the vorticity transport equation are solved in a 
coupled manner, there is no need to identify which boundary 
condition is used for which equation. 


A Least Squares Formulation 
Consider the functional 


* (v-q-s) + | vxq - w | 2 dV 



Minimizing I with respect to q yields 


5 1 » 0 = JJJ (v-q - s) v»5q + (vxq - Z) • vx«q dV 
V 

- -JJJ v(7.q - s) • sq + vx(vxq - Z) • 5q dV 
V 

+ JJ (v.q - s) Sq-n + (vxq - Z) • (nxsq) dA 
A 


A more detailed analysis of this equation is shown in the Appendix. 

Thus, the Euler-Lagrange equation associated with minimizing I is 
identical to equation (8), provided the boundary terms vanish. The 
latter condition is satisfied if i) the velocity vector q is specified at 
the boundary; ii) the tangential velocity components are specified 
and equation (1) is imposed at the boundary; or iii) the normal 
velocity component is specified and the tangential components of 
equation (2) are imposed at the boundary. 

Assume the boundary conditions of iii) are used to determine the 
velocity field. To solve the vorticity transport equations, three 
boundary conditions are required; these are the two tangential 
velocity components prescribed at the boundary as well as the 
compatibility condition on the vorticity field . 

It is noticed by Gunzberger and Patterson (10) that the latter is a 
natural boundary condition associated with the Galerkin formulation 
of the vorticity transport equation. Moreover, taking the divergence 
of the vorticity equation, it is shown that v uJ=o is governed by a 
linear homogeneous equation and hence v-u3 vanishes everywhere. 

To avoid the coupling between the velocity components at the 
curved boundaries, the boundary condition i) may be used instead, 
to determine the velocity field. The boundary conditions for the 
vorticity transport equations are still 

uxn * vxqxn. and v*o> - 0 

Special treatment may be required, however, to determine accurately 
the tangential vorticity component at a curved boundary in terms of 
the derivatives of the velocity components. 



Numerical Techniques and Discr etization 


Finite Difference Discretization. 

The terms are discretized using central difference schemes and 
staggered grid formulation. This formulation defines vorticity to be 
in the center of the grid while the velocity components are located 
along the edges. The extension of this approach to 3-D schemes is 
also applied. The method of solution is as follows: For the 2-D 

formulation, the equations are solved fully coupled using Newton's 
method of linearization and a direct solver package (UNPACK) . An 
outline of what subroutines are called and their roles are given in the 
Appendix. For the 3-D formulation, the Poisson equations are solved 
directly. The LU decomposition for the velocity components are only 
done once and can be stored for later use. This will help with the cpu 
requirement. The vorticity transport equations are solved using a 
zebra line relaxation method for each of the 3 planes. The 2-D 
staggered grid arrangement is show in Fig. 1 while the 3-D 
arrangement is show in Fig. 2. 


Finite Element Discretization. 

A discrete approximation of the functional I can be easily 
constructed using finite element techniques. Using standard bilinear 
shape functions for the velocity components, the functional I can be 
evaluated in terms of the nodal values. Upon minimization over 
each element and assembling the contribution from all elements, the 
nodal equations are readily obtained. 

A Galerkin method due to Swartz and Wendroff (11) ( see also 
Fletcher (12) ) is used for the vorticity equation. The boundary 
conditions are implemented as discussed before. There resulting 
nonlinear system of algebraic equations are solved by Newtons 
method and a direct solver. 

For a simple geometry, a finite difference discretization over a 
staggered grid leads to the five point stencil for the Poisson 
equations. The stencil for the finite element method is shown in Fig. 
1 . 



Finite Volume Discretization ( Generalized Coordinate System ). 

When arbitrary geometries are considered the choice of the 
discretization scheme is crucial. The discretized governing equations 
and their compatibility relations must be satisfied as well a certain 
geometrical identities when a generalized coordinate system is 
employed. The finite-volume approach can give accurate 
conservative approximations as was pointed out by Vinokur (13). 
This method was also appled by Rosenfeld et al (14) using a 
primitive variable formulation on a staggered mesh and a 
generalized non-orthogonal system using the velocity/vorticity 

formulation. 

The choice of the dependent variables is also critical in obtaining 
the accurate solutions. The scaled contravariant velocity 
components in a staggered grid arrangement are used as the 
dependent variables instead of the conventional Cartesian velocity 
components. This choice is essential for mass conservation on the 
discrete level. 


The Integral Formulation of the Governing equations are: 
jf q .n dS * 0 ( mass conservation ) 

}} nxqdS = JJJ cadV (definition of vorticity) 
fj <*n dS » 0 ( conservation of vorticity) 

\\ ( vq n ♦ n x u ) dS = 0 

[} (nq«-n«q-n»)dS = 0 


The geometric identity 
ff ndS - C 

imposes the condition that the cell is closed. Another conditon that 
must be also satisfied is that the sum of the cell volumes must equal 
the total 

volume of the flow region. 

These conditions are automatically satisfied for the special case of 
staggered Cartesian coordinates for a cubic driven cavity. 


Numerical Results 

The flow over a backward facing step is simulated using exactly 
the same geometry as Ref. 2. The domain of integration is bounded 
by solid surfaces at the top and the bottom of the backward facing 
step. At the inlet station, a parabolic velocity profile is imposed 
between the shoulder and the upper surface and it is assumed that 
the flow is parallel( i.e. the upstream effect of the shoulder is 
neglected). At the exit, the flow is assumed to be fully developed 
(independent of x) ,and so v=0. The tangential velocity component 
(u) is not known and the continuity equation is enforced at the exit 
boundary. Fig. 3 shows the problem formulation for this case. 

Fig. 4 shows the computing time for the Cray X-MP as compared to 
the INS3D code and the fractional step method. The streamlines are 
plotted in Figures 5-8 for the Reynolds number ranges of 133 
through 1867. 


Experimental resuts not only gave the expected primary zone of 
recirculation region attached to the step but also showed additional 
recirculation region downstream of the steps and on both sides of the 
channel. Experiments have shown that transition to turbulence 
happens at around RE=900. This means that the laminar and steady 
flow assumption that was mentioned earlier is no longer valid. 
However, much higher Reynolds number are analyzed to show that 
the method is robust. Fig. 1 1 shows the convergence history for both 


finite difference and finite element formulation. Quadratic 
convergence is obtained independent of the grid size and Re provided 
the initial guess is available. Comparison with other numerical 
results and experimental data are also given. Mesh refinement is 
also performed and it is seen that the finite element results are less 
sensitive to the grid size. 

Fig. 12 shows the convergence history for the 3-D driven cavity 
problem. 

The domain is bounded by solid surfaces everywhere with the top 
wall moving at a uniform velocity of U=1 in the positive x - direction. 

Grid sizes of 9x9x9 up to 17x17x17 are used for Re=100. These 
conditions are chosen to obtain direct comparison with the published 
results. Figure 13 shows the distribution of the u-velocity 
component along the centerline(z=0.5,y=0.5,x=0.5) as compared with 
numerical results obtained by Dennis et. al. 

Figure 14 shows the variation of the velocity component as a 
function of y at z=0. and x=0.5 as compared with Dennis et al. 

The vorticity contours compare well with Osswald et al for the 
same grid size of 17x17x17. The contours for station 
x=0.5,y=0.5,z=0.5 are shown in Fig 15 while the contours for station 
x=0.7 1785, y=0.71785,z=0.7 1785 are show in Fig. 16. 

The current work is on 2-D unsteady flow around a cylinder using 
curvilinear coordinates. 
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Discretization 


Finite difference: -discretize by central difference 

-staggered grid 



w=vorticity 
u,v=vel comp 


Finite Element : -discretize using bilinear shape 

functions 

"regular grid 



Uj(x,y) = a Uj + b uj x + c uj y *d u . xy 


Figure 1 





Proposed Method for 3-D 
Velocitv-Vorticitv Formulation 


Governing equations (3-d) 

(11) U XX + Uyy + U ZZ = (W2) Z - (Wg)y 

< 12 > v xx + v yy + V 2Z = ( w 3 )x ‘ ( W 1 )y 

(13) W xx + Wyy+ W 2Z = (W-J )y- (W 2 ) * 

(14) U(w^) x +V(Wi )y + w(w-| ) Z -W-| U x -W2Uy-WgU z = 

1_[ W 1xx + w 1yy+ w 1zz ] 

Re 

(15) u(w2) x +v(W2)y +w(w2) z -w-| v x -W2Vy-W3V z = 

1[ W 2xx +W 2 yy +W 2zz ] 

Re 

(16) u(w 3 ) x +v(w 3 )y +w(w 3 ) z -w-| W x -W 2 Wy-W 3 W z 

l [w 3xx +w 3yy +w 3zz] 

Re 


(17) W 1= Wy-v z 

(18) w 3 = U 2' W X 

(19) w 3 = v x -Uy 


3-D staggered grid arrangement: 


y-dlr 



cell 1 jk 
— z-dir 


v-rtir 


Figure 2. 



Problem Formulation (2-D1 


(8) 

V 2 U = - Wy 

( 9 ) 

V^v = w x 

(10) 

(UW) X + (VW)y 


w 


Re 


yy 


non-dimensionalized variables : 

x = x'/h , y = y'/h , u = u’/Umax , v = v'/Umax 

w=w'/(Umax/h) , Re = (2/31Umax(2h1 

v 

Umax = maximum inlet velocity 
h = step height 
L = channel length 
v = kinematic viscosity 

Test model (2-D backward facing step ) : 



Figure 3 



COMPUTING TIME FOR BACKWARD FACING STEP 
ON CRAY X-MP 


( T= t cpu sec/grid point /Iteration) 


Re 

INS3D 

FRACTIONAL STEP 

PRESENT METHOD 

100. 

2.0x10-6 

1.5x10-6 

*cpu lime range of 

200. 

2.0x10-6 

3.1 x 1 0-6 

1.5x10-6 to 3.0x10-6 

300. 

2.6x10-6 

5.2 xlO-6 

for Re=100 to 1800 

400. 

6.7x10-6 

8.2x10-6 


500. 

1.0 xlO-5 

1.1 xlO-5 


600. 

1.0x10-5 

1.2x10-5 


700. 

.97 x 10-5 

1.3 x 1 0-5 


800. 

10 x 10-5 

1.4x10-5 



INS3D and Fractional Step data obtd. from S.Rogers,D.KwaM.Chang, 
Numerical Solution of Incompressible NS Eo ns In 3-D General 
Curvilinear Coord, ( NASA TM 86840 ). Jan 1986 


Figure 4. 



staggered grid formulation 



STREAMLINE CONTOURS 105x31 


Re=133. 


finite element formulation 



STREAMLINE CONTOURS 105x31 


FigureB 


staggered grid formulation 



STREAMLINE CONTOURS 105x31 


Re= 600. 


finite element formulation 



S T R E A 


MLINE CONTOURS 105x31 


Figure fr 



■ taggored grid formulation 



STREAMLINE CONTOURS 105x31 


Re=1200. 


finite ole m e n t for m u 1 a t i o n 



0.0 10.0 20.0 30.0 40.0 

STREAMLINE CONTOURS 105x31 


Figure 7, 




Re=l867. 


finite element formulation 



STREAMLINE CONTOURS 105x31 


Figure 8 
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Convergence history of 
a finite element calculation 


Convergence history of 
a finite difference calculation 
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RE = 100. 3-D DRIVEN CAVITY (17X17X17) GRID 
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rl= u-vel 
r2= v-vel 
r3=w-vel 


r4= wl transport eqn 
r5= w2 transport eqn 
r6= w3 transport eqn 


r6= dlv( q ) 
r7= dlv( w ) 
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Figure 12. Convergence History 
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RE=100. UftRIflTION UF U AT Z = fl.. V=.Q 
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Figure 16 Contours of Normal Velocity Components at (x= 0.7 1 075, 
y=0.7 1 075, z=0.71875 ) 




